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Abstract 

Coexistence properties of the hard-core attractive Yukawa potential with inverse-range parameter 
K = 9, 10, 12 and 15 are calculated by applying canonical Monte Carlo simulation. As previously 
shown for longer ranges, we show that also for the ranges considered here the coexistence curves 
scaled by the critical density and temperature obey the law of corresponding states, and that a 
linear relationship between the critical density and the reciprocal of the critical temperature holds. 
The simulation results are compared with the predictions of the self-consistent Ornstein-Zernike 
approximation, and a good agreement is found for both the critical points and the coexistence 
curves, although some slight discrepancies are present. 
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I. INTRODUCTION 



Short-range potential models of simple analytic form have been the object of intense in- 
vestigation over the last years. The main reason for such an interest is that they provide an 
approximate representation of the effective interactions and phase behavior experimentally 
observed in a variety of real systems. For instance, globular protein solutions, colloidal sus- 
pensions, and fullerenes at high temperature 

Q-B- This scenario has been further enriched 
by recent investigations of gas-solid, solid-solid, and a metastable vapor-liquid transition 
in one of the most studied short-range models, namely the hard-core attractive Yukawa 
(HCAY) potential ^ 



12|. The HCAY model is given by 



u{r) 



(1) 



oo, if r < cr, 

— ecrexp[— /t(r — cr)]/r, if o" < r, 
where e is the depth of the potential and a is the hard-sphere diameter. The range of the 
potential is tuned by varying k. Our results are given in dimensionless units, such that 
r* = r/a for distance, T* = ksT/e for temperature [ks is the Boltzmann's constant), 
p* = pcr^ for density. 

Many theoretical approaches have been used for studying the HCAY fluid, for example: 
the Barker-Henderson perturbation theory [2], [l3|, the thermodynamic perturbation theory 



(TPT) 0, Q-17| and a variety of other mean field approaches based on truncation or 
approximate summation of perturbative series Q, |l9], the mean spherical approximation 
(MSA) HI 20l422(|. the modified hypernetted chain approximation (MHNC) [23|] and other 

25! . the density functional theory js]. 



integral equation theories 



24 



261. the self-consistent 



Ornstein-Zernike approximation (SCOZA) jll. I271 3Cll| . and the hierarchical reference theory 
(HRT) 23|, l3l|] , among others. For all theories, the most difficult regime to deal with is that 
of short interaction range, where an accurate localization of the vapor-liquid coexistence 
curve and of the critical point presents severe difficulties. Sometimes, the shape of the 
coexistence curve is not reproduced correctly and critical points are not located correctly 
either, especially the critical density. We found that the critical density data reported by 
these approaches show different tendencies for shorter-range attractions, when we plot p* 
versus 1/T*. These issues were some of the motivations to carry out this work. 

For large values of k > 7, the computer simulation data for HCAY fluid are quite 
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scarce. In this regime, Gibbs ensemble Monte Carlo (GEMC) cannot reproduce binodal 
cu.ve= BQ. The pha.e diagra,n at large . wa= =tud.ed by Monte Carlo supplemented by 

thermodynamic integration (MC-TI) in Ref. [7|, but the results do not include the critical 
density. In our previous works we reported the coexistence properties of HCAY fluid for long 
range tails. These properties were compared to those reported in the literature, and some 
differences were found among them In a recent work, Singh [36| has reported the 



coexistence curves of short-range attractive Yukawa (SR-HCAY) fluid with k = 8 — 10, using 
Grand-canonical transition-matrix Monte Carlo (GC-TMMC) with the histogram reweight- 
ing method. However, this method proves difficult to use in order to calculate coexistence 
properties at low temperature and high liquid densities [36|. So, NVT-MC simulations ap- 
pear to be a more efficient method to calculate such properties for very short-range systems, 
where these conditions are met. 

In view of the above considerations, the main goals of this paper are to demonstrate that 
there is a linear dependence between the critical density and the reciprocal of the critical 
temperature even for short interaction ranges, where theoretical approaches show different 
tendencies, and to report a systematic study of the liquid-vapor phase diagrams of SR- 
HCAY fluid with K = 9, 10, 12 and 15, using canonical ensemble Monte Carlo simulation, 
where most theoretical approaches fail and other simulation techniques cannot reproduce 
the binodal curves. Finally, once again we have conflrmed that the coexistence curves of SR- 
HCAY model follow the law of corresponding states, i.e., curves corresponding to different 
K fall on top of each other within a high degree of accuracy, provided the density p and 
temperature T are rescaled by their critical values pc, Tc. 

Besides, we have also reported the binodal curves and critical points predicted by SCOZA 
for the above values of n. Although the application of SCOZA to narrow Yukawa potentials 
has been considered before [l|, |23|, an extensive comparison with simulation data for the 
phase diagram in this regime has not been possible so far because of the aforementioned 
paucity of simulations at large k, and the present study represents a good opportunity to 
perform it. Results were obtained both by the standard version of SCOZA considered in 



Refs. 



28|, and by a modifled version developed in 29|], where consistency with the 



virial route, which is disregarded in the original implementation, is partially taken into 
account. The comparison shows that some discrepancies between SCOZA and simulations 
are deflnitely present: speciflcally, the standard formulation of SCOZA overestimates both 
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the critical density and critical temperature, while the converse happens with the modified 
formulation. The latter also predicts a nearly constant critical density at values of k, for 
which the simulation show that this quantity still changes significantly. Nevertheless, the 
overall agreement with the simulations is quite satisfactory. 

The paper is organized as follows: in Sec. Hi] the details of our simulations are briefiy 
summarized; in Sec. Illlt an overview of both the standard and the modified versions of 
SCOZA is given, and some issues relevant for the application of the theory to short-range 
interactions are discussed; in Sec. |IV]our results are presented and discussed; in Sec. IVlour 
conclusions are drawn. 



II. SIMULATION DETAILS 

Canonical ensemble Monte Carlo (NVT-MC) simulations of HCAY fiuid interface have 
been performed using Metropolis algorithm. To prepare the initial configuration we placed 

= 1726 particles in a face-centered cubic array in the middle of the simulation cell, which 
allowed us to obtain two interfaces with a vapor phase surrounding the liquid when the 



system is equilibrated (see Fig 2 of Ref . [33| ) 



Simulations were performed on a parallelepiped cell with dimensions = Ly = 10, and 
Lz = 40. Additional MC runs were carried out with = Ly = 14 and a larger number of 
to check for finite-size effects, and no such effect was found. The cutoff 



particles 
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radius was selected to Rc = 2. Periodic boundary conditions were applied in all three 
directions and the neighbor list was applied to speed up calculation. The maximum particle 
displacement was adjusted to give a 45% acceptance rate. The simulations were performed in 
cycles; in every cycle it was attempted to move all the particles that are in the neighbor list. 
The system was equilibrated for 10^ cycles and the coexistence properties were calculated for 
10^ cycles divided into 100 blocks of 10^ cycles each one, in order to avoid density fluctuation. 

The density proflles, p{z), were calculated every 50 cycles and at the end of each sim- 
ulation run, every proflle was fltted to a hyperbolic tangent function to obtain py and p2 



values 



33 



35| . The critical parameters for the HCAY fluid were calculated by using the 



rectilinear diameter law 39|] and the universal value Pcoex = 0.325 for the critical exponent 
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Pcoex which gives the curvature of the hquid-vapor coexistence curve in the critical region. 
All the critical point parameters we calculated are listed in Table HI 

III. THEORY 

A detailed description of SCOZA and its application to a Yukawa fluid has been given 
elsewhere and we refer the interested reader to the previous literature on this subject 



23 
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28| . Here we recall that, as usual in integral-equation theories, in the SCOZA 
the Ornstein-Zernike (OZ) equation is closed by an approximate ansatz which involves the 
direct correlation function c(r) and the radial distribution function g{r). In the standard 
formulation of SCOZA, one requires that i) the core condition g{r) = Q be satisfied inside 
the hard core r < a; and ii) the contribution to c(r) outside the hard core due to tail of the 
potential be linear in the potential itself. For a HCAY potential one has then: 

qir) = if r < 0" , 

c(r) = CHsir) + -ft'exp[— K(r — cr)]/r if cr < r , 

where CHs{f) is the direct correlation function of the hard-sphere fluid. The functional 
form of this closure is similar to that of well-established theories such as the mean spherical 
approximation (MSA), where CHs{f^) is identically vanishing for r > cr, or the optimized 



random-phase approximation (ORPA) j40|. However, in these approaches the amplitude 
K of the direct correlation function coincides with the inverse temperature /3 = l/{kBT), 
while in the SCOZA is an a priori unknown function of the thermodynamic state, to be 
determined so that consistency between the compressibility and the internal energy route to 
the thermodynamic is obtained. The consistency requirement is embodied in the following 
condition on the reduced compressibility Xred and the excess internal energy per unit volume 
u: 

where Xred and u are obtained via the compressibility and internal energy route respectively. 
Eq. (|3]) together with closure ([2]) give a closed partial differential equation (PDE) for m, 
that can be solved numerically, provided the initial condition at /3 = and the boundary 
conditions at low and high density are specified. The solution procedure is made easier by 
the fact that, for a given the OZ equation with closure (|2]) can be solved analytically, 
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if c(r) for r > o" is given by a superposition of Yukawa tails [4l|]. This is indeed the case, 
provided one adopts the Waisman parametrization for the hard-sphere direct correlation 
function CHs{f) [42]. According to this, Cnsif^) for r > a is also given by a Yukawa tail 
with known, density-dependent amplitude and range, so that c(r) for r > a consists of the 
superposition of two Yukawas. On the other hand, for non- Yukawa intermolecular potentials 
the SCOZA c(r) for a < r is not of Yukawa form, irrespective of the parametrization adopted 
for CHs{f)i and in order to implement SCOZA a fully numerical solution of the OZ equation 



is usually required [43|, . 



In conventional SCOZA, the virial route is not included in the consistency condition. In 
order to enforce consistency among all the three routes, one has to endow Eq. ([2]) with one 



more degree of freedom. To this purpose, the following form for c(r) was proposed 29l |: 



c(r) = CHs{r) + /3exp[-K(r - cr)]/r 

-|- if exp[— 2;(r — cr)]/r if a < r , (4) 

where the amplitude of the (attractive) Yukawa potential has now been set to (5 as in the 
MSA, while both the amplitude H and the inverse range z of the extra Yukawa tail must 
be determined by imposing consistency among virial, compressibility, and internal energy. 
This approach is very similar to that originally proposed in Ref. [45!], which is obtained 
from Eq. (jlj) by having the last Yukawa to account also for the hard-sphere contribution 
to c(r) outside the repulsive core, thereby setting cusir) = in Eq. (jlj). A fully numerical 



implementation of that closure was considered in Ref. 23|, but the numerical algorithm failed 
to converge in the critical region. In a recent work 29|, a different strategy was adopted: on 
the one hand, the analytical solution of the OZ equation for a c(r) given by a superposition 
of Yukawas (three when also c//5(r) is taken into account) for a < r was again exploited. 
On the other hand, consistency with the virial route was required to hold only at the critical 
point, and z was fixed at the value at which this condition is satisfied, thereby making the 
theory considerably easier to implement. This is the approach used here together with the 
standard SCOZA of Eqs. (E]), ([3]). In the following the standard and modified versions of 
SCOZA will be referred to as S-SCOZA and M-SCOZA respectively. 

Before concluding this Section, it is worthwhile adding some remarks about the relevance 
of the high-density boundary condition for the SCOZA PDE. While the boundary condition 
at low density is trivial, since one must have Xred{p = ^i P) = 1, w(p = 0,/3) = for every /3, 



the behavior at high density poses more problems. In fact, this is not known beforehand, 
so that there is a certain amount of arbitrariness on the form of the corresponding bound- 
ary condition. For instance, in Refs. 28|, the boundary condition was imposed on 



d'^u/dp'^, which was set to the value obtained by the so-called high-temperature approxima- 
tion (HTA) |40j, i.e., by replacing u with its zeroth-order contribution in an expansion in 
powers of (3. In Ref. 29|, instead, the boundary condition was imposed directly on u rather 
than on its derivative, and u was determined using the g{r) obtained via the ORPA, which 
amounts to setting K = 13 in Eq. (|2]). Since there is not any obvious recipe for singling out 
a specific high-density boundary condition, a natural requirement is that the results should 
not depend on the detailed form of such a condition, or on the density po at which it is im- 
posed. This turns out to be indeed the case, provided po is high enough. For Lennard- Jones 
like tails, such as that corresponding to the widely adopted value 2; = 1.8 setting the 
high-density boundary at Pq = 1 proves to be amply sufficient. However, for interactions of 
shorter range, such as those considered in Ref. [l| as well as in the present paper, po must be 
shifted to substantially higher densities. This was also observed in a fully numerical imple- 
mentation of SCOZA for a square- well fluid 43||. In fact, for the largest value of z considered 
here, z = 25 (see Table [T]), po must be at least as high as the close-packing density pg = \/2, 
and at even larger z one needs to move po beyond close packing. This requirement sounds 
quite unphysical, but one should keep in mind that here the effect of the singular repulsion is 
described at the level of the Carnahan-Starling equation of state 40|], which is itself unaware 
of close packing, and has the compressibility vanish at the unphysical value p* = 6/n rather 
than at p* = V2. The HTA for u is expected to become exact at close packing [46], but in 
view of the above, it is not surprising that one might need to go beyond close packing for 
the hard-sphere repulsion to overwhelm the attractive part of the interaction in g{r). 

If the boundary is not located at sufficiently high po, the results turn out to be very 
sensitive with respect to both a shift in po, and the details of the boundary condition. This 
observation is relevant for some of the results reported in Ref. 29|], where both S-SCOZA 
and M-SCOZA calculations were performed for different values of z, but the high-density 
boundary was always kept fixed at pg = 1. This is adequate for 2; < 8, but for larger z it leads 
to an underestimation of both the critical density and the critical temperature with respect 
to the genuine SCOZA result, which becomes more and more severe as z is increased. The 
results presented here are not affected by this spurious effect, and are meant to supersede 
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those shown in Table I of Ref. 



29| 



IV. RESULTS AND DISCUSSIONS 



We start our discussion with an analysis of the critical values of HCAY fluid. Table [T] 
presents the reduced critical temperature T* and density p* for a wide range of k predicted by 
both simulations and several theoretical approaches, namely, the S-SCOZA and M-SCOZA 



described in the previous section, the TPT by Zhou and the equation of state (EOS) by 
Duh and Mier-Y-Teran [22]. We recall that, as A is increased, the liquid- vapor critical point 
becomes metastable with respect to freezing. According to many theoretical and simulation 



studies 



for A > 7 the metastable regime has alreday set in. 



K 


^ C 




Method 


Ref. 


1.8 


1.180 


0.313 


NVT-MC 


m 


1.8 


1.180 


0.315 


GC-TMMC 


[36j 


1.8 


1.212 


0.312 


GC-FSS" 


[28j 


1.8 


1.219 


0.314 


S-SCOZA 


[28j 


1.8 


1.246 


0.310 


TPT 


M 


1.8 


1.240 


0.318 


MSA-EOS 


[22j 


2.0 


1.050 


0.322 


NVT-MC 


m 


2.0 


1.088 


0.323 


S-SCOZA 


this work 


2.5 


0.840 


0.336 


NVT-MC 


m 


2.5 


0.871 


0.342 


S-SCOZA 


this work 


3.0 


0.721 


0.356 


NVT-MC 


m 


3.0 


0.722 


0.355 


GC-TMMC 


[36j 


3.0 


0.740 


0.359 


S-SCOZA 


this work 


3.0 


0.764 


0.379 


MSA-EOS 


[22j 


4.0 


0.581 


0.380 


NVT-MC 


[33j 


4.0 


0.572 


0.385 


GC-TMMC 


[36J 


4.0 


0.591 


0.389 


S-SCOZA 


[23j 


4.0 


0.579 


0.375 


M-SCOZA 


this work 


4.0 


0.588 


0.414 


TPT 


M 
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4.0 


0, 


.614 


0.428 


MSA-EOS 


fool 

[22] 


5.0 


0, 


.500 


0.39 


NVT-MC 


fool 

[33] 


5.0 


0, 


.509 


0.416 


S-SCOZA 


this work'' 


5.0 


0, 


.497 


0.393 


M-SCOZA 


this work 


5.0 


0, 


.530 


0.472 


MSA-EOS 


rr»f-»l 

[22j 


6.0 


0, 


.448 


0.412 


NVT-MC 


To nl 

[33J 


6.0 


0, 


.456 


0.438 


S-SCOZA 


this work'' 


6.0 


0, 


.445 


0.407 


M-SCOZA 


this work 


7.0 


0, 


.414 


0.422 


NVT-MC 


fool 

[33] 


7.0 


0, 


.419 


0.457 


S-SCOZA 


[23] 


7.0 


0, 


.409 


0.418 


M-SCOZA 


this work 


7.0 


0, 


.411 


0.502 


TPT 


[14] 


8.0 


0, 


.382 


0.447 


GC-TMMC 


[36] 


8.0 


0, 


.392 


0.474 


S-SCOZA 


this work 


8.0 


0, 


.383 


0.426 


M-SCOZA 


this work 


9.0 


0, 


.365 


0.456 


ATT IV /r/~i 

NVT-MC 


this work 


9.0 


0, 


.362 


0.454 


GC-TMMC 


To *-il 

[36] 


9.0 


0, 


.370 


0.489 


S-SCOZA 


this work 


9.0 


0, 


.368 


0.433 


M-SCOZA 


this work 


10.0 


0, 


.347 


0.466 


ATT nn TV /r/~l 

NVT-MC 


this work 


10.0 


0, 


.343 


0.471 


GC-TMMC 


To 

\M 


10.0 


0, 


.352 


0.503 


S-SCOZA 


this work 


10.0 


0, 


.345 


0.439 


TV T O / \ r7 A 

M-SCOZA 


this work 


10.0 


0, 


.361 


0.652 


TV (TO A T — ■, /' \ ( ~\ 

MSA-EOS 


[22] 


11.0 


0, 


.337 


0.515 


S-SCOZA 


this work 


11.0 


0, 


.331 


0.444 


M-SCOZA 


this work 


12.0 


0, 


.323 


0.480 


NVT-MC 


this work 


12.0 


0, 


.324 


0.525 


S-SCOZA 


this work 


12.0 


0, 


.319 


0.447 


M-SCOZA 


this work 


15.0 


0, 


.292 


0.515 


NVT-MC 


this work 



15.0 


0.293 


0.551 


S-SCOZA 


this work 


15.0 


0.291 


0.455 


M-SCOZA 


this work 


25.0 


0.235 


0.555 


MC-TI 




25.0 


0.236 


0.600 


S-SCOZA 


this work 


25.0 


0.241 


0.459 


M-SCOZA 


this work 


25.0 


0.239 


0.632 


TPT 


M 



TABLE I: Critical parameters of HCAY Fluid. 
'^Gran canonical Monte Carlo with finite-size scaling. 
''The calculation for this value of k had already been per- 
formed in Ref. [l], but the result had not been tabulated. 



'^Pc was obtained using Eq. (4) from Ref. 



35|. 



The critical density as a function of the inverse of the critical temperature as obtained 
from our simulations and the aforementioned theories is shown in Fig. [TJ The critical simu- 
lation data obtained by Singh using GC-TMMC, also displayed in the figure, are in excellent 
agreement with ours, which shows that both techniques give similar results within the error 
bar. At small k, all theories give similar results, and are in good agreement with the simu- 
lations. In fact, in the long-range regime even a simple mean-field approach along the lines 
of the van der Waals theory provides an accurate description of the system, which has been 
shown to become exact in the limit where the range of the tail potential goes to infinity, and 



its strength goes to zero 47|]. As the interaction range gets shorter, discrepancies become 



more and more pronounced, especially as far as p* is concerned. In particular, the data for 



p* reported by Duh and Mier-Y-Teran [22| using their equation of state (EOS) for HCAY 
fiuid strongly overestimate the simulations. The TPT by Zhou jl^ and the S-SCOZA show 
a better agreement with the simulation data. The latter are bracketed by the S-SCOZA 
and M-SCOZA results: S-SCOZA overestimates both p* and T*, while the converse is true 
for M-SCOZA. For both approaches, the relative error with respect to the simulations is in 
general considerably smaller for T* than for p*. 

The most remarkable feature of the simulation results is the linear trend for p* vs. 1/T* 
up to inverse-range parameters as large as k = 15, corresponding to 1/T* = 3.425. A similar 
behavior is displayed by S-SCOZA, as well as by TPT. As k is further increased, the slopes 
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0.675 



0.575 - 



0.475 - 



0.375 - 



ONVT-MC 
+ GC-TMMC 
□ S-SC0ZA(1.9) 
VM-SC0ZA(1.9) 

TPT 
>EOS 




0.5 



1.5 



2.5 

1/T„ 



3.5 



4.5 



FIG. 1: Relation between critical density and inverse critical temperature of HCAY model for 
different simulation and theoretical results. NVT-MC with 1.8 < k < 15 from this and previous 
Sa-S^l- GC-TMMC with 1.8 < k < 10 [36]. SCOZA with 1.8 < k < 25 from this and 



works 



i],y,i28i 



previous works [11,123,123], M-SCOZA with 4 < k < 25 from this work, EOS with 1.8 < k < 12 
and TPT with 1.8 < k < 25 Lf 



of both SCOZA and TPT plots start to decrease, and deviations from linearity become 
apparent. This is consistent with the fact that, as the potential range goes to zero and 1/T* 
diverges, one expects p* to tend to a finite, non-vanishing value. Such a qualitative behavior 
is predicted by Baxter's analytical solution of his adhesive hard-sphere (AHS) model 



and has recently been supported by MC simulations of the Ah 



S model itself 



49 



48|, 
5oi as 



well as of square- well fluids with very short attraction range [5l|- It may also be recalled 

ue of p* in the limit of vanishing attraction range is 



that for the HCAY fluid, a finite va^ 
found analytically also within MSA 52|]. We did not perform simulations for n > 15, but 
we took the critical temperature for k = 25 obtained in Ref . [3] , and inserted it into Eq. (4) 



of Ref. 



351 ] to obtain the critical density. The result, also plotted in Fig. HI shows a trend 



similar to that of S-SCOZA. In this respect, it may be observed that, although M-SCOZA 
has the best agreement with simulations for k, < 15, it overemphasizes the tendency to 



saturation at larger n 



531 ] ■ which is better described by S-SCOZA. 



We now turn to the dependence of the critical temperature on the interaction range. 
For hard-core plus attractive tail potentials, it has been suggested that the second virial 
coefficient at the critical temperature i?2(T'c) has a practically constant value i?2, and that 
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coefficient B2(T) is 



54 



55|. Specifically, 



short-range tails follow a universal behavior, provided the second viria 
used as the temperature-like variable instead of the temperature itself 
it should be possible to map the system into a hard-core plus square-well fluid with the 
same hard-core diameter, and a temperature-dependent range 6 determined so as to give 
the same B2 at the same reduced temperature T* as those of the original interaction. As 
a consequence, for small enough ranges the reduced critical temperature should fulffil the 
relation 



IP* 



1 + 



(5) 



1 + 5)3-1^ 

where B2^ is the virial coefficient of the hard-sphere fluid. Subsequent work [sl, [ssl [sgI [s^I 
has shown that, even for short-range interactions, i?2(Tc) changes systematically with the 
attraction range. Nevertheless, at short range Eq. ([5]) does represent quite accurately our 
results for T*. This is shown in Fig. |2l where T* as predicted by both simulations and 
theories has been plotted as a function of the effective square-well range 6, and compared 
with Eq. (j5]) by taking B2 as a ffiting parameter. It can be noticed that all the theories 
considered here agree reasonably well with the simulations, and the discrepancies are smaller 
than those found for p*. Moreover, for 6 < 0.2, T* is indeed well represented by Eq. ([5]) 
with B^/B^^ ~ -1.35. In fact, at short range the effective range 6 is much less sensitive 
than B2(T) to small changes in temperature, so that Eq. (|5]) can be satisfied to a very good 
degree, even though i?2(Tc) is not actually constant. 

If B2{Tc) tends to a finite limit B^"^ for 5 ^ 0, Eq. with B^ = B^"^ will hold for 
potentials of arbitrarily short range down to the AHS limit. Such a behavior was found by 
Largo et al. [51] for the square- well fiuid. We are not able to state if this is the case also for 
the Yukawa fiuid since, for the interval of k considered in this work, our simulation results 
do not provide a clear evidence of B2{Tc) saturating to a finite limit. In this respect, we 
should notice that the value -Bj/-^!^"^ — —1-35 found here is lower than both the result of 



Largo et al. [51] B^/B^^ = -1.174 and that of Miller and Frenkel |49| B^/B^^ = -1.21. 

Our new phase equilibrium data of HCAY fiuid with k = 9, 10, 12, and 15 are shown 
in Fig. [31 and compared to those reported by Singh [36| for k = 9 and n = 10. An 
excellent agreement is found between them. Therefore, both simulation techniques are good 
for calculating such properties. It is worth noting that the NVT-MC technique is more 
efficient for the calculation of the coexistence properties of very short-range systems and 
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FIG. 2: Relation between critical temperature and effective range 6 (see text) of the HCAY model 
for different simulation and theoretical results. Symbols have the same meaning as in Fig. [TJ MC- 

n 

TI with 3.9 < K < 100 from Ref. [7]. The red line represents the prediction of Eq. ([5]). The inset 
is an enlargement of the short-range region. 
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FIG. 3: Vapor-liquid coexistence curves of SR-HCAY fluid, 
this work. Other symbols are GC-TMMC data from Ref. 



Open symbols are NVT-MC data from 



36|. Dotted lines are S-SCOZA results 



from this work. Starbursts and circles are S-SCOZA and M-SCOZA critical points, respectively. 

at low temperature, where other simulation techniques such as GEMC or GC-TMMC run 
into considerable difficulties, because the high liquid density of the systems makes particle 
insertion very time-consuming. Our simulation results for the coexistence curve are reported 
in Table HIl Figure E] also shows the coexistence curves and the critical points predicted 
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K 




PL 


Pv 


9.0 


0.330 


0.945 


0.0650 




0.335 


0.915 


0.0766 




0.340 


0.890 


0.0900 




0.345 


0.850 


0.1100 




0.350 


0.810 


0.1510 




0.355 


0.741 


0.2110 


10.0 


0.315 


0.965 


0.0621 




0.320 


0.930 


0.0834 




0.325 


0.895 


0.1005 




0.330 


0.860 


0.1233 




0.335 


0.820 


0.1638 


12.0 


0.300 


0.945 


0.1021 




0.305 


0.905 


0.1282 




0.308 


0.873 


0.1454 




0.310 


0.855 


0.1603 




0.312 


0.832 


0.1802 




0.315 


0.802 


0.2022 


15.0 


0.280 


0.960 


0.1201 




0.282 


0.928 


0.1380 




0.284 


0.901 


0.1550 




0.286 


0.855 


0.1900 




0.288 


0.805 


0.2310 



TABLE II: Coexistence properties of SR-HCAY fluid at different interaction ranges. 

by S-SCOZA, and the M-SCOZA critical points. M-SCOZA coexistence curves have not 
been reported, because we could not rule out the possibility that they may be affected by 
significant numerical errors. It appears that the main source of discrepancy between the 
simulation and theoretical phase diagrams lies in the position of the critical points whose 
densities and temperatures, as observed above, are both slightly overestimated by S-SCOZA. 
Figure m shows the reduced density as a function of the reduced temperature. All curves of 
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FIG. 4: Reduced vapor-liquid coexistence curves of SR-HCAY fluid. NVT-MC data from this work 



and Ref. BJ. GC-TMMC data from Ref. 



short range HCAY fluid form a single master curve within a very good degree of accuracy, 
as it was reported by one of us in a previous work for longer range [sS]. As the range 
is decreased, a slight trend towards a widening of the coexistence curve can be detected, 
but this effect remains small for the ranges considered here. Also, the coexistence curves 
reported by Singh 36| using GC-TMMC are on the master curve. The same behavior has 
been found for the Mie(n,m) potential 58N60|. These results confirm that, even for the 
short range interactions considered here, with inverse-range parameter as large as k = 15, 
the vapor-liquid phase diagrams of HCAY fluid obey the law of corresponding states. This 
fact could be used to test theoretical approaches. 



V. CONCLUSIONS 



We have carried out a systematic study of the liquid-vapor phase diagrams for hard- 
core attractive Yukava potential with k, = 9, 10, 12 and 15, using NVT-MC simulation. 
Our coexistence curves for k = 9 and 10 were compared to those reported by Singh using 
GC-TMMC. An excellent agreement was found between them. Besides, we have shown 
that reduced coexistence curves of HCAY fluid from k = 7 to k = 15 follow the law of 
corresponding states, as it was shown in a previous work [35j for longer interaction ranges. 
Once again, we have confirmed the linear relationship between the critical density and inverse 
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critical temperature for critical simulation data in this interval of k. The simulation results 
were compared with several theories, including the standard formulation of SCOZA (S- 
SCOZA) as well as a modified version (M-SCOZA) implementing a partial consistency with 
the virial route, which is disregarded in the standard version. The comparison shows that 
S-SCOZA slightly overestimates both the critical temperature and the critical density with 
respect to the simulation results, while the converse is true for M-SCOZA. Nevertheless, the 
agreement with the simulations is quite good for all the values of k considered in this paper. 
Both simulation and theory give a dependence of the critical temperature as a function of the 
interaction range in substantial agreement with Noro-Frenkel scaling. However, whether the 
second virial coefficent at the critical temperature B2{Tc) tends to a constant as k diverges, 
thereby making Noro-Frenkel scaling hold at arbitrarily short range as in the square- well 
case 5l|], could not be established on the basis of our simulation data, and remains to be 
assessed. 
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